pacman::p_load(rsample)


boot_percentile_function <- function(alpha, 
                                     boots, 
                                     statistics) {
  alpha <- alpha
  statistics <- statistics
  boots <- boots
  bootstrap_percentiles <- purrr::map(
    .x = boots,
    .f = ~ rsample::int_pctl(
      .data = .x,
      alpha = alpha, 
      statistics = statistics
    )
  )
}


# Define Variables  ####

set.seed(853147)

outcome_vars <- c("news_sharing", 
                  "news_distinguish")
experimental_vars <- c("fastpaced_nudge", 
                       "accuracy_nudge")


# Bootstrap of Fastpaced x Accuracy x Authoritarian Models  ####

tic()
## bootstraps
bootstraps <- data_nested %>%
  ungroup() %>% 
  select(data_type, outcome, data) %>% 
  # subset the data for computational feasibility
  mutate(data = map(data,
                    ~ select(.,
                             any_of(outcome_vars),
                             any_of(authoritarian_vars), 
                             any_of(experimental_vars), 
                             any_of(observational_vars),
                             starts_with("num_"), 
                             matches("index"), 
                             "country_id", 
                             "ID"))) %>%
  # nest data by respondent ID for clustered bootstrap resampling
  mutate(data = map(data, 
                    ~ nest(.x, data = -c("ID")))) %>% 
  #sample 5000 times over each subset of the data
  mutate(boots = map2(
    .x = data,
    .y = outcome, 
    .f = ~ rsample::bootstraps(.x,
                               strata = .x[['country_id']],
                               times = boot_times,
                               apparent = T
    )
  )) %>% 
  # estimate model on each of the splits
  mutate(boots = map2(
    .x = boots,
    .y = outcome, 
    ~ mutate(
      .x,
      # lm fastpaced accuracy authoritarian interact simple
      lm_fastpaced_accuracy_authoritarian_simple_tidy = pmap(
        list(..1 = .x$splits, 
             ..2 = .y), 
        .f = ~ map(authoritarian_vars,
                   function(authoritarian_vars) {
                     fit_lm_fastpaced_accuracy_authoritarian(
                       outcome = ..2, 
                       data = analysis(..1) %>% 
                         unnest(data),  
                       model = "simple", 
                       tidy = TRUE, 
                       authoritarian_measure = paste(authoritarian_vars))
                     }
                   )
        ),
      # lm fastpaced accuracy authoritarian interact full
      lm_fastpaced_accuracy_authoritarian_full_tidy = pmap(
        list(..1 = .x$splits, 
             ..2 = .y), 
        .f = ~ map(authoritarian_vars,
                   function(authoritarian_vars) {
                     fit_lm_fastpaced_accuracy_authoritarian(
                       outcome = ..2, 
                       data = analysis(..1) %>% 
                         unnest(data),  
                       model = "full", 
                       tidy = TRUE, 
                       authoritarian_measure = paste(authoritarian_vars))
                   }
                )
        )
      )
    )) %>% 
  select(data_type, outcome, boots) %>% 
  mutate(boots = map(boots, 
                     ~ select(., c("splits", "id", matches("tidy"))))
  )


# Confidence Intervals of Fastpaced x Accuracy x Authoritarian Models  ####

cis_fastpaced_accuracy_authoritarian <- bootstraps %>%
  select(boots, data_type, outcome) %>%   
  # create a column that stores the wide tibbles of the nested list of tibbles
  # preserves the rset structure in boots which is necessary for pctl_int
  mutate(lm_fastpaced_accuracy_authoritarian_simple = map(boots, 
                                                 ~ unnest_wider(., 
                                                                col = lm_fastpaced_accuracy_authoritarian_simple_tidy) %>% 
                                                   select(-splits, -id, -matches("lm_")) %>% 
                                                   rename_with(.cols = everything(), 
                                                               ~ paste0(.x, "_simple")))) %>% 
  # add the wide tibbles back to the rset object (boots)
  mutate(boots = map2(.x = boots, 
                      .y = lm_fastpaced_accuracy_authoritarian_simple, 
                      .f = ~ bind_cols(.x, .y) %>% 
                        # deselect nested list of tibbles
                        select(., -lm_fastpaced_accuracy_authoritarian_simple_tidy))) %>% 
  # deselect helper column with wide tibbles
  select(-lm_fastpaced_accuracy_authoritarian_simple) %>%   
  # create a column that stores the wide tibbles of the nested list of tibbles
  # preserves the rset structure in boots which is necessary for pctl_int
  mutate(lm_fastpaced_accuracy_authoritarian_full = map(boots, 
                                               ~ unnest_wider(., 
                                                              col = lm_fastpaced_accuracy_authoritarian_full_tidy) %>% 
                                                 select(-splits, -id, -matches("lm_")) %>%
                                                 rename_with(.cols = everything(), 
                                                             ~ paste0(.x, "_full")))) %>% 
  # add the wide tibbles back to the rset object (boots)
  mutate(boots = map2(.x = boots, 
                      .y = lm_fastpaced_accuracy_authoritarian_full, 
                      .f = ~ bind_cols(.x, .y) %>% 
                        # deselect nested list of tibbles
                        select(., -lm_fastpaced_accuracy_authoritarian_full_tidy))) %>% 
  # deselect helper column with wide tibbles
  select(-lm_fastpaced_accuracy_authoritarian_full) %>% 
  # 90 CI 
  mutate(boot_pctl_90_simple =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.10,  
              statistics = paste0(authoritarian_vars, "_simple"))
          }
    ))) %>%   
  # 95 CI
  mutate(boot_pctl_95_simple =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.05,  
              statistics = paste0(authoritarian_vars, "_simple"))
          }
    ))) %>% 
  # 99 CI 
  mutate(boot_pctl_99_simple =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.01,  
              statistics = paste0(authoritarian_vars, "_simple"))
          }
    ))) %>% 
  # 90 CI 
  mutate(boot_pctl_90_full =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.10,  
              statistics = paste0(authoritarian_vars, "_full"))
          }
    ))) %>%   
  # 95 CI
  mutate(boot_pctl_95_full =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.05,  
              statistics = paste0(authoritarian_vars, "_full"))
          }
    ))) %>% 
  # 99 CI 
  mutate(boot_pctl_99_full =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.01,  
              statistics = paste0(authoritarian_vars, "_full"))
          }
    ))) %>% 
  select(data_type, outcome, starts_with("boot_pctl")) %>% 
  pivot_longer(cols = starts_with("boot_pctl"), 
               names_to = c("ci_level", "model_specification"), 
               names_pattern = "boot_pctl_(.*)_(.*)",
               values_to = "boot_pctl") %>% 
    unnest(boot_pctl) %>% 
    mutate(authoritarian_specification = names(boot_pctl)) %>% 
    unnest(boot_pctl)


cis_fastpaced_accuracy_authoritarian %>% 
  write_csv(.,
            file = "data/bootstrapped_cis/cis_fastpaced_accuracy_authoritarian.csv")

toc()


# By Country: Bootstrap of Fastpaced x Accuracy x Authoritarian Models  ####
rm(bootstraps, 
   cis_fastpaced_accuracy_authoritarian)



tic()
## bootstraps
bootstraps <- data_nested_by_country %>%
  ungroup() %>% 
  select(country, data_type, outcome, data) %>% 
  # subset the data for computational feasibility
  mutate(data = map(data,
                    ~ select(.,
                             any_of(outcome_vars),
                             any_of(authoritarian_vars), 
                             any_of(experimental_vars), 
                             any_of(observational_vars),
                             starts_with("num_"),
                             matches("index"), 
                             "country_id", 
                             "ID"))) %>%
  # nest data by respondent ID for clustered bootstrap resampling
  mutate(data = map(data, 
                    ~ nest(.x, data = -c("ID")))) %>% 
  #sample 5000 times over each subset of the data
  mutate(boots = map2(
    .x = data,
    .y = outcome, 
    .f = ~ rsample::bootstraps(.x,
                               times = boot_times,
                               apparent = T
                               )
  )) %>% 
  # estimate model on each of the splits
  mutate(boots = map2(
    .x = boots,
    .y = outcome, 
    ~ mutate(
      .x,
      # lm fastpaced authoritarian interact simple
      lm_fastpaced_accuracy_authoritarian_simple_tidy = pmap(
        list(..1 = .x$splits, 
             ..2 = .y), 
        .f = ~ map(authoritarian_vars,
                   function(authoritarian_vars) {
                     fit_lm_fastpaced_accuracy_authoritarian(
                       outcome = ..2, 
                       data = analysis(..1) %>% 
                         unnest(data),  
                       model = "simple", 
                       tidy = TRUE, 
                       authoritarian_measure = paste(authoritarian_vars), 
                       by_country = TRUE)
                   }
        )
      ),
      # lm fastpaced accuracy authoritarian interact full
      lm_fastpaced_accuracy_authoritarian_full_tidy = pmap(
        list(..1 = .x$splits, 
             ..2 = .y), 
        .f = ~ map(authoritarian_vars,
                   function(authoritarian_vars) {
                     fit_lm_fastpaced_accuracy_authoritarian(
                       outcome = ..2, 
                       data = analysis(..1) %>% 
                         unnest(data),  
                       model = "full", 
                       tidy = TRUE, 
                       authoritarian_measure = paste(authoritarian_vars), 
                       by_country = TRUE)
                   }
        )
      )
    )
  )) %>% 
  select(country, data_type, outcome, boots) %>% 
  mutate(boots = map(boots, 
                     ~ select(., c("splits", "id", matches("tidy"))))
  )




# By Country: Confidence Intervals of fastpaced x Authoritarian Models ####

cis_fastpaced_accuracy_authoritarian_by_country <- bootstraps %>%
  select(country, boots, data_type, outcome) %>%   
  # keep only fastpaced interaction models
  mutate(boots = map(.x = boots, 
                     .f = ~ select(., splits, id, matches("authoritarian")))) %>%   
  # create a column that stores the wide tibbles of the nested list of tibbles
  # preserves the rset structure in boots which is necessary for pctl_int
  mutate(lm_fastpaced_accuracy_authoritarian_simple = map(boots, 
                                                 ~ unnest_wider(., 
                                                                col = lm_fastpaced_accuracy_authoritarian_simple_tidy) %>% 
                                                   select(-splits, -id, -matches("lm_")) %>% 
                                                   rename_with(.cols = everything(), 
                                                               ~ paste0(.x, "_simple")))) %>% 
  # add the wide tibbles back to the rset object (boots)
  mutate(boots = map2(.x = boots, 
                      .y = lm_fastpaced_accuracy_authoritarian_simple, 
                      .f = ~ bind_cols(.x, .y) %>% 
                        # deselect nested list of tibbles
                        select(., -lm_fastpaced_accuracy_authoritarian_simple_tidy))) %>% 
  # deselect helper column with wide tibbles
  select(-lm_fastpaced_accuracy_authoritarian_simple) %>%   
  # create a column that stores the wide tibbles of the nested list of tibbles
  # preserves the rset structure in boots which is necessary for pctl_int
  mutate(lm_fastpaced_accuracy_authoritarian_full = map(boots, 
                                               ~ unnest_wider(., 
                                                              col = lm_fastpaced_accuracy_authoritarian_full_tidy) %>% 
                                                 select(-splits, -id, -matches("lm_")) %>%
                                                 rename_with(.cols = everything(), 
                                                             ~ paste0(.x, "_full")))) %>% 
  # add the wide tibbles back to the rset object (boots)
  mutate(boots = map2(.x = boots, 
                      .y = lm_fastpaced_accuracy_authoritarian_full, 
                      .f = ~ bind_cols(.x, .y) %>% 
                        # deselect nested list of tibbles
                        select(., -lm_fastpaced_accuracy_authoritarian_full_tidy))) %>% 
  # deselect helper column with wide tibbles
  select(-lm_fastpaced_accuracy_authoritarian_full) %>% 
  # 90 CI 
  mutate(boot_pctl_90_simple =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.10,  
              statistics = paste0(authoritarian_vars, "_simple"))
          }
    ))) %>%   
  # 95 CI
  mutate(boot_pctl_95_simple =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.05,  
              statistics = paste0(authoritarian_vars, "_simple"))
          }
    ))) %>% 
  # 99 CI 
  mutate(boot_pctl_99_simple =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.01,  
              statistics = paste0(authoritarian_vars, "_simple"))
          }
    ))) %>% 
  # 90 CI 
  mutate(boot_pctl_90_full =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.10,  
              statistics = paste0(authoritarian_vars, "_full"))
          }
    ))) %>%   
  # 95 CI
  mutate(boot_pctl_95_full =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.05,  
              statistics = paste0(authoritarian_vars, "_full"))
          }
    ))) %>% 
  # 99 CI 
  mutate(boot_pctl_99_full =  map(
    .x = boots, 
    ~ map(authoritarian_vars, 
          function(authoritarian_vars) {
            rsample::int_pctl(
              .data = .x,
              alpha = 0.01,  
              statistics = paste0(authoritarian_vars, "_full"))
          }
    ))) %>% 
  select(country, data_type, outcome, starts_with("boot_pctl")) %>% 
  pivot_longer(cols = starts_with("boot_pctl"), 
               names_to = c("ci_level", "model_specification"), 
               names_pattern = "boot_pctl_(.*)_(.*)",
               values_to = "boot_pctl") %>% 
  unnest(boot_pctl) %>% 
  mutate(authoritarian_specification = names(boot_pctl)) %>% 
  unnest(boot_pctl)


cis_fastpaced_accuracy_authoritarian_by_country %>% 
  write_csv(.,
            file = "data/bootstrapped_cis/cis_fastpaced_accuracy_authoritarian_by_country.csv")


rm(bootstraps, 
   cis_fastpaced_accuracy_authoritarian_by_country)

toc()